library(CausalGAM)


M <- 1000

n.vec <- c(250, 500, 1000)

output.mat.a <- matrix(NA, M*3*3*2*7, 8)
output.mat.i <- matrix(NA, M*3*3*2*7, 6)
output.mat.r <- matrix(NA, M*3*3*2*7, 6)

counter <- 1
for (n in n.vec){
  for (alpha.ind in 1:3){
    for (beta.ind in 1:2){
      for (i in 1:M){
        filename <- paste("../../data/MCdata.n", n, ".alpha", alpha.ind,
                          ".beta", beta.ind,
                          ".rep", i, ".Rdata", sep="")
        
        load(filename)

        for (spec in 1:7){
          cat("n = ", n, "\n")
          cat("alpha = ", alpha.ind, "\n")
          cat("beta = ", beta.ind, "\n")
          cat("i = ", i, "\n")
          cat("spec = ", spec, "\n\n")
          
          if (spec == 1){
            out <- estimate.ATE(pscore.formula=x~lo(z1,z2),
                                outcome.formula.t=y~lo(z2)+lo(z3),
                                outcome.formula.c=y~lo(z2)+lo(z3),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1,
                                data=mydata, divby0.tol=0, nboot=0)            
          }
          else if (spec == 2){
            out <- estimate.ATE(pscore.formula=x~lo(z1,z2)+lo(z3),
                                outcome.formula.t=y~lo(z1)+lo(z2)+lo(z3),
                                outcome.formula.c=y~lo(z1)+lo(z2)+lo(z3),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1,
                                data=mydata, divby0.tol=0, nboot=0)
          }
          else if (spec == 3){
            out <- estimate.ATE(pscore.formula=x~lo(z1),
                                outcome.formula.t=y~lo(z2),
                                outcome.formula.c=y~lo(z2),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1.0,
                                data=mydata, divby0.tol=0, nboot=0)
          }
          else if (spec == 4){
            out <- estimate.ATE(pscore.formula=x~lo(z2),
                                outcome.formula.t=y~lo(z3),
                                outcome.formula.c=y~lo(z3),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1.0,
                                data=mydata, divby0.tol=0, nboot=0)
          }
          else if (spec == 5){
            out <- estimate.ATE(pscore.formula=x~lo(z2),
                                outcome.formula.t=y~lo(z2),
                                outcome.formula.c=y~lo(z2),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1.0,
                                data=mydata, divby0.tol=0, nboot=0)
          }
          else if (spec == 6){
            out <- estimate.ATE(pscore.formula=x~lo(z2)+lo(z3),
                                outcome.formula.t=y~lo(z1, z2),
                                outcome.formula.c=y~lo(z1, z2),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1.0,
                                data=mydata, divby0.tol=0, nboot=0)
          }
          else if (spec == 7){
            out <- estimate.ATE(pscore.formula=x~lo(z2),
                                outcome.formula.t=y~lo(z2)+lo(z3),
                                outcome.formula.c=y~lo(z2)+lo(z3),
                                pscore.family=binomial(probit),
                                treatment.var="x",
                                var.gam.plot=FALSE,
                                outcome.family=gaussian,
                                variance.smooth.span=-1.0,
                                data=mydata, divby0.tol=0, nboot=0)
          }

          
          output.mat.a[counter,1] <- out$ATE.AIPW.hat
          output.mat.a[counter,2] <- n
          output.mat.a[counter,3] <- alpha.ind
          output.mat.a[counter,4] <- i
          output.mat.a[counter,5] <- spec
          output.mat.a[counter,6] <- beta.ind
          output.mat.a[counter,7] <- out$ATE.AIPW.sand.SE
          output.mat.a[counter,8] <- (abs(out$ATE.AIPW.hat - 5.00) / out$ATE.AIPW.sand.SE) < qnorm(.975)

          output.mat.i[counter,1] <- out$ATE.IPW.hat
          output.mat.i[counter,2] <- n
          output.mat.i[counter,3] <- alpha.ind
          output.mat.i[counter,4] <- i
          output.mat.i[counter,5] <- spec
          output.mat.i[counter,6] <- beta.ind

          output.mat.r[counter,1] <- out$ATE.reg.hat
          output.mat.r[counter,2] <- n
          output.mat.r[counter,3] <- alpha.ind
          output.mat.r[counter,4] <- i
          output.mat.r[counter,5] <- spec
          output.mat.r[counter,6] <- beta.ind

          
          
          counter <- counter + 1
        } ## end spec loop


        
      }
    }
  }
}

colnames(output.mat.a) <- c("ATE.hat",
                            "n", "alpha.ind", "rep",
                            "spec", "beta.ind",
                            "ATE.sand.se", "cover.sand")

colnames(output.mat.i) <- c("ATE.hat",
                            "n", "alpha.ind", "rep",
                            "spec", "beta.ind")

colnames(output.mat.r) <- c("ATE.hat",
                            "n", "alpha.ind", "rep",
                            "spec", "beta.ind")


write.csv(output.mat.a, file="../../results/AIPW.csv")
write.csv(output.mat.i, file="../../results/IPW.csv")
write.csv(output.mat.r, file="../../results/Reg.csv")
